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Abstract 

Using data from the FOCUS (E831) experiment at Fermilab, we present a new 
measurement of the branching ratio for the Cabibbo-favored decay mode D° 
K-K-K+TT+. 

From a sample of 143±19 fully reconstructed K~ K~ K'^tt'^ events, wc mea- 

sure r{D^ K-K-K+Tr+)/T{D^ if-7r-7r+7r+) = 0.00257 ± 0.00034(siat.) ± 
0.00024(syst.)- 

A coherent amplitude analysis has been performed to determine the resonant 
substructure of this decay mode. This analysis reveals a dominant contribution 
from (f> and K (892) states. 



1. Introduction 

Hadronic decays of charm mesons have been extensively studied in recent 
years. Dalitz plot analyses of D meson decays in three-body final states have 
revealed a rich resonant substructure, showing a dominance of quasi two- 
body modes. However, much less information is available on the resonant 
substructure of decays with more than three final state particles. These multi- 
body D decays account for a large fraction of the hadronic decay width. If 
the amplitude analyses of multi-body decays confirm the picture drawn by 
the Dalitz plot analyses, then one can make a more complete comparison 
with theoretical models, which have been developed mainly to describe the 
two-body and quasi-two-body decay modes [1,2,3,4,5]. 

We present a new study of the D° — > K~K''K'^7r'^ decay using data from 
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the FOCUS experiment. This is an interesting decay mode: although Cabibbo 
favored, it is strongly suppressed by phase-space and requires the production 
of an ss pair, either from the vacuum or via final state interactions (FSI). 

We measure the branching ratio ^r{D^K-Tr-^+'K+) perform a coherent 

amplitude analysis to determine its resonant substructure. 

FOCUS is a charm photoproduction experiment [6] which collected data dur- 
ing the 1996-97 fixed-target run at Fermilab. Electron and positron beams 
(with typically 300 GeV endpoint energy) obtained from the 800 GeV/c Teva- 
tron proton beam produce, by means of bremsstrahlung, a photon beam which 
interacts with a segmented BeO target. The mean photon energy for triggered 
events is ~ 180 GeV. A system of three multicell threshold Cerenkov coun- 
ters performs the charged particle identification, separating kaons from pions 
with momenta up to 60 GeV/c. Two systems of silicon microvertex detectors 
are used to track particles: the first system consists of 4 planes of microstrips 
interleaved with the experimental target [7]; the second system consists of 12 
planes of microstrips located downstream of the target. These detectors pro- 
vide high resolution in the transverse plane (approximately 9 /xm), allowing 
the identification and separation of charm primary (production) and secondary 
(decay) vertices. The charged particle momenta are determined by measuring 
their deflections in two magnets of opposite polarity through flve stations of 
multiwire proportional chambers. 



2. Analysis of the decay mode ^ K K K'^n'^ 

The final states are selected using a candidate driven vertex algorithm [6]. A 
secondary vertex is formed from the four candidate tracks. The momentum 
of the resultant candidate is used as a seed track to intersect the other 
reconstructed tracks to search for a primary vertex. The confidence levels of 
both vertices are required to be greater than 1%. Once the production and 
decay vertices are determined, the distance i between them and its error a£ 
are computed. The quantity £/cr£ is an unbiased measure of the significance 
of detachment between the primary and secondary vertices. This is the most 
important variable for separating charm events from non-charm, prompt back- 
grounds. Signal quality is further enhanced by cutting on Iso2, which is the 
confidence level that other tracks in the event might be associated with the 
secondary vertex. To minimize systematic errors on the measurements of the 
branching ratio, we use identical vertex cuts on the signal and normalizing 
mode, namely i/a£ > 6, and Iso2 < 1 %. We also require the D° momentum 
to be in the range 25-250 GeV/c (a very loose cut) and the primary vertex to 
be formed with at least two reconstructed tracks in addition to the D° seed. 

The only difference in the selection criteria between the two decay modes lies 
in the particle identification cuts. The Cerenkov identification cuts used in 
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FOCUS are based on likelihood ratios between the various particle identifi- 
cation hypotheses. These likelihoods are computed for a given track from the 
observed firing response (on or off) of all cells within the track's (/3 = 1) 
Cerenkov cone for each of our three Cerenkov counters. The product of all 
firing probabihties for all the cells within the three Cerenkov cones produces 
a x^-like variable Wi = — 2 In(Likelihood) where i ranges over the electron, 
pion, kaon and proton hypotheses [8]. All kaon tracks are required to have 
Ak = — Wk (kaonicity) greater than 2, whereas the pion tracks are re- 
quired to have A^r = Wk — W^tt (pionicity) exceeding 0.5. 



Using the set of selection cuts just described, we obtain the invariant mass 
distributions for K~ K~ K^ti~ and K~ti~ti^ti^ shown in Fig. 1. In Fig. la 
the K~ K~ K~^TT^ mass plot is fit with two Gaussians with the same mean but 
different sigmas to take into account the change in resolution with momentum 
of our spectrometer [6] plus a second-order polynomial. A log-hkelihood fit 
returns a signal of 143 ± 19 D° — > K~K~K'^7r~ events. The large statistics 
K~7r~TT~^TT~^ mass plot of Fig. lb is fitted in the same way (two Gaussians 
plus a second-order polynomial). The fit gives a signal of 64576 ± 360 
K^Tf^TT^n^ events. 



3. Relative Branching Ratio 



The evaluation of relative branching ratios requires yields from the fits to 
be corrected for detection efficiencies, which differ among the various decay 
modes due to differences in both spectrometer acceptance (due to different 
Q values and resonant substructure of the two decay modes) and Cerenkov 
identification efficiency. 



From Monte Carlo simulations, we compute the relative efficiency to be: 

e{D"^K-K-K+TT+) 

r{D^^K-K-K+7T+) 



e{D"^K K K+n+) _ q ^ Q.OIO. Combined with the yields measured above 



we find a branching ratio of wr.o~^ ^- _ + \ J = 0.00257 ± 0.00034. 

Our final measurement has been tested by modifying each of the vertex and 
Cerenkov cuts individually. The branching ratio is stable versus several sets of 
cuts as shown in Fig. 2 (we vary the confidence level of the secondary vertex 
from 1% to 10%, Iso2 from 10"^ to 1, i/at from 5 to 15, from 0.5 to 4.5 
and finally A^ from 0.5 to 4.5). 



Systematic uncertainties on branching ratio measurements can come from dif- 
ferent sources. We determine four independent contributions to the systematic 
uncertainty: the split sample component, the fit variant component, the com- 
ponent due to the particular choice of the vertex and Cerenkov cuts (discussed 
previously), and the limited statistics of the Monte Carlo. 
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The split sample component takes into account the systematics introduced 
by a residual difference between data and Monte Carlo, due to either a pos- 
sible mismatch in the reproduction of the momentum or the change in 
the experimental conditions of the spectrometer during data collection. This 
component has been determined by splitting data into four independent sub- 
samples, according to the momentum range (high and low momentum) and 
the configuration of the vertex detector, that is, before and after the insertion 
of an upstream silicon system[7]. A technique employed in FOCUS and in 
the predecessor experiment E687, modeled after the S-factor method from the 
Particle Data Group [9], is used to try to separate true systematic variations 
from statistical fluctuations. The branching ratio is evaluated for each of the 
4 statistically independent sub-samples and a scaled variance is calculated a 
(the errors are boosted when /{N — 1) > 1). The split sample variance (Jsput 
is defined as the difference between the reported statistical variance and the 
scaled variance, if the scaled variance exceeds the statistical variance [10]. 

Another possible source of systematic uncertainty is the fit variant. This com- 
ponent is computed by varying, in a reasonable manner, the fitting conditions 
on the whole data set. In our study we fixed the widths of the Gaussians to 
the values obtained by the Monte Carlo simulation, changed the background 
parametrization (varying the degree of the polynomial), and used one Gaus- 
sian instead of two. In addition we considered the variation of the computed 
efficiency, both for — > K~ K~ K~^7r~ and the normalizing decay mode, due 
to the different resonant substructure simulated in the Monte Carlo. The BR 
values obtained by these variants are all a priori likely; therefore this uncer- 
tainty can be estimated by the rms of the measurements [10]. 

Analogous to the fit variant, the cut component is estimated using the stan- 
dard deviation of the values obtained from the many sets of cuts shown in 
Fig. 2. Actually this is an overestimate of this component because the event 
samples of the various cut sets are different. 

Finally, there is a contribution due to the hmited statistics of the Monte Carlo 
simulation used to determine the efficiencies. The resulting systematic errors 
are summarized in Table 1. Adding in quadrature the four components, we 
obtain the total systematic error also shown in Table 1. 

The final branching ratio result is shown in Table 2 along with a comparison 
to previous measurements. 



4. Amplitude analysis of ^ K K i^+7r+ 

A fully coherent amplitude analysis was performed to determine the reso- 
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Source 


Percent 


hpiit sample 


U.U/o 


Fit Variant 


/I 1 07 
4.1/0 


Set of cuts 


8.3% 


MC statistics 


1.2% 


Total systematic 


9.4% 



Table 1 

Contribution in percent to the systematic uncertainties of the branching ratio 
r{D^ K-K-K+'K+)/T{D^ ii'-7r-7r+7r+). 



Experiment 


V{D'^^K-K-K+-K+) 


Events 


T(D^-^K--K--K+'K+) 


E687 [11] 
E791 [12] 

FOCUS (this result) 


0.0028 ± 0.0007 ± 0.0001 
0.0054 ± 0.0016 ± 0.0008 
0.00257 ± 0.00034 ± 0.00024 


20 ±5 
18.4 ±5.3 

143 ± 19 



Table 2 

Branching ratio measurement and comparison with other experiments. 

nant substructure of the — > K~ K~ K~^7r~^ decay. While a large number 
of intermediate states could lead to the K~K~K'^7r'^ final state, phase space 
hmitations restrict the possible contributions. A plot of the K~K^ invariant 
mass shows a clear contribution (Fig. 3, top left plot). It is also possible, 
in principle, to have K~ contribution via /o(980) and ao(980). There are 
large uncertainties in the line shape of these resonances and in their coupling 
to K~K^. They do not appear to be required by the fit and are therefore not 
included. 

Contributions from ^(800)-^'"-^'"'" and K {lAd>0)K^ might also be present. 
However, these resonances are broad scalar states, with no characteristic an- 
gular distribution that could distinguish them from the non-resonant mode, 
with the present level of statistics. 

Modes containing a K (892) can also contribute, even though the nominal 
K (892) mass is just above the kinematical limit of the K~t{^ spectrum. Since 
K (892) is a narrow vector meson, its contribution can be distinguished by 
the angular distribution of the decay products, even without a clear mass 
peak. 

We consider four decay amplitudes: non-resonant, (j)K (892), (j)K~TT^ and 
K {892)K-K+. In our default model, model A, we use all four amplitudes. 
We also consider a simpler model, model B, which does not include the two 
K (892) amphtudes. 
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The formalism used in this amphtude analysis is a straightforward extension 
to four-body decays of the usual Dalitz plot fit technique. The is a spin 
zero particle, so the four-body decay kinematics are defined by five degrees of 
freedom. 

Individual amplitudes, Aj., for each resonant mode are constructed as a prod- 
uct of form factors, relativistic Breit-Wigner functions, and spin amplitudes 
which account for angular momentum conservation. We use the Blatt-Weisskopf 
damping factors [14], Fi, as form factors {I is the orbital angular momentum 
of the decay vertex). For the spin amplitudes we use the Lorentz invariant 
amplitudes [13], which depend both on the spin of the resonance(s) and the 
orbital angular momentum. The relativistic Breit-Wigner is 



m? — rriQ -\- imQT{m) ' 



where 

^ \PoJ ^sQ 



In the above equations m is the two-body invariant mass, mo and s are the 
resonance nominal mass and spin, and p* — p*{m) is the breakup momentum 
at resonance mass m. 

Since there are two identical kaons, each amplitude is Bose-symmetrized. 
The overall signal amplitude is a coherent sum of the individual amplitudes, 

A = J2kCkAk, assuming a constant complex amplitude for the non-resonant 
mode. The coefficients are complex numbers to be determined by the ffi. 
The overall signal amplitude is corrected on an event-by-event basis for the 
acceptance, which is nearly constant across the phase space. 

In the amplitude analysis we have taken events having a K~ K~ K^ti^ in- 
variant mass within the interval Md ± 10 MeV/c^. In this interval there are 
139 signal and 65 background events. The finite detector resolution causes a 
smearing of the edges of the five-dimensional phase space. This effect is ac- 
counted for by multiplying the overall signal distribution by a Gaussian factor, 
g{M), where M is the K~ K~^K~7r~^ mass. The normalized signal probability 
distribution is 

P5(0) = ^e(0)p(0)^(M) I J2c,AM |^ 

where represents the coordinates of an event in the five dimensional phase 
space, £(0) is the acceptance function, and is the phase space density. 
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Two types of background events were considered: random 0's combined with 
a K^Ti^ pair, and random combinations of K~ K^t^^ . Inspection of the 
side bands of the K~K~K~^'k^ mass spectrum indicate that nearly 30% of the 
background events are of the former type. We assume the random K~K~K'^7r'^ 
combinations to be uniformly distributed in phase space, while for the (p back- 
ground we assume an incoherent sum of Breit-Wigners with no form factors 
and no angular distribution. The overall background distribution, kept fixed 
in the fit, is a weighted, incoherent sum of these two components. The overall 
background distribution is also corrected for the acceptance (assumed to be 
the same as for the signal events) on an event-by-event basis, and multiplied 
by an exponential function b{M), to account for the detector resolution. The 
normalized background probability distribution is 

An unbinned maximum likelihood fit was performed, minimizing the quantity 
w = — 21n£. The hkelihood function, £, is 

n [Psm+PBm 

events 

Neither the acceptance function nor the phase space density depend on the fit 
parameters c^, so the term — 2/n [e(0)p(0)] is irrelevant to the minimization. 
The acceptance correction is important only for the normalization integrals 
Ns and Nb- 

Decay fractions are obtained from the coefficients determined by the fit, 
after integrating the overall signal amplitude over the phase space: 



/*; — 





1 CkAk p 


I 


Id<f>\ 


J2j CjAj 


2 



We fit the data to both models A and B. We find that model B, with only the 
non-resonant and (f)K~Tr~^ contributions, does not provide a good description 
of the data. The inchision of the two K (892) amplitudes results in a much 
better fit, with an improvement in Aw of 138. 

Results from the fit with model A are shown in Table 3. The dominant con- 
tribution comes from (f)K (892). Adding the contributions from the 
two K (892) amplitudes, we see that they account for nearly 70% of the total 
decay width. This is somewhat surprising, given the very small phase space. In 
Fig. 3 the K^K~ and K^n^ projections of events used in the amplitude anal- 
ysis are superimposed on the fit result (top two plots). The projections from 
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the background model are shown in the shaded histograms. The remaining 
plots are two-dimensional projections of the events in the signal region. 

Wc have performed a log-likelihood test to check our ability to distinguish 
between models A and B. Two ensembles of 10,000 mini-MC samples were 
generated, one simulated according to model A and another according to 
model B. For each sample in each ensemble we compute the quantity Aw — 
2 ln>C^ — 2 luvCs, where jCb and jCa are the likehhoods calculated with models 
B and A, respectively. The resulting Aw distributions are shown in Fig. 4. On 
the right we see the At;; distribution computed from the model A ensemble, 
and, on the left. Aw with the model B ensemble. The two distributions are well 
separated showing that we can easily distinguish between these two models. 
Moreover, the value of Aw obtained from the real data (138) is consistent with 
the distribution from model A, showing that this is indeed a better description 
of the data. 

The goodness-of-fit was assessed in two ways. We have estimated the confi- 
dence level of the fit using a test. The five invariants used to define the 
kinematics of this decay are the four K~K^ and K~ti:^ masses squared, plus 
either the K~K~ or K'^n'^ mass squared. Due to the limited statistics we 
have integrated over the latter invariant and divided the other four into two 
bins, yielding a total of sixteen cells. A was computed and the estimated 
confidence level was 35%. The confidence level obtained with model B was 
6x10"^^ 

Given the limited statistics we have also estimated the confidence level using 
a method which is often less stringent than the x^- this second method the 
confidence level is estimated using the distribution oi w — —2 In £ from an 
ensemble of mini-MC samples generated with the parameters of model A. This 
distribution is approximately Gaussian. A confidence level can be estimated 
by the fraction of samples in which the value of w exceeds that of the data. 
With this technique we estimate a GL of 86% for the fit with model A. 

As in the branching ratio measurement, we consider split sample and fit variant 
systematic uncertainties, using for the former the same sub-samples described 
previously. The dominant contributions from fit variant systematic errors 
come from variations on the signal amplitudes (removing the Blatt-Weisskopf 
form factors and replacing the non- resonant amplitude by K {14:S0)K'^K~), 
variations in the background relative fractions, and variations of the analy- 
sis cuts. Split sample and fit variant errors are added in quadrature. Table 3 
shows the Model A fit results. The systematic errors are the second errors in 
the fractions and phases. 
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mode 


Tnf^.ffTiit.iirlf 


nhaspf dp£f 1 


fraction (%^ 




1 
1 


n 

u 


ZIZ U ZiZ 1 


(l)K~'K+ 


0.60 ± 0.12 


194 ± 24 ± 8 


18 ± 6 ± 4 


K*\892)K+K- 


0.65 ± 0.13 


255 ± 15 ± 4 


20 ± 7 ± 2 


non-resonant 


0.55 ± 0.14 


278 ± 16 ± 42 


15 ± 6 ± 2 



Table 3 

Results from the best fit (Model A). The second error on the fractions and phases 
is systematic. 

5. Conclusions 

Using data from the FOCUS (E831) experiment at Fermilab, we have studied 
the Cabibbo-favored decay mode K~K~K~^7r~^. 

A comparison with the two previous determinations of the relative branching 

ratio ^Y{D°^K-n-n+n+) ^^^^^ imprcssivc improvement in the accuracy of 
this measurement. 

A coherent amphtude analysis of the K~K~K'^7r^ final state was performed 
for the first time, showing that the dominant contribution comes from D° — > 
4>K (892). This channel, D° decaying to two vector mesons, corresponds to 

50% of the L)0 K-K-K+Tr+ decay rate. The Z*°(892) amplitudes, in 
spite of the limited phase space, account for nearly 70% of the total decay 
width. Looking at the K^K~ spectrum we see that over 60% comes from 
r>o ^ 0X*°(892) and (j)K- 7r+. 
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Fig. 1. Invariant mass distribution for K~K~K^TT~^{a) and /'r^7r^7r+7r+(b). The fit 
(solid curve) is to two Gaussians for the signal plus a quadratic polynomial for the 
background. 
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Fig. 2. Branching ratio T{D° K' R- K"^ -k^) /T {D^ iC-7r-7r+7r+) versus sev- 
eral sets of cuts. Prom left to right, we vary the confidence level of the secondary- 
vertex from 1% to 10% (10 points), Iso2 from 10~^ to 1 (7 points), L / (Jl from 5 
to 15 (11 points), Ak from 0.5 to 4.5 (9 points) and A^r from 0.5 to 4.5 (9 points). 
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Fig. 3. The top two plots are projections of the invariant mass distribution K~ 
and K'-K^ , with the fit results superimposed (solid histograms). The shaded his- 
tograms are the background projections. The left scatter plots show the data pop- 
ulation for K~K^ and K~tt~^ from the signal region, while the right ones show the 
fit results. Due to two identical particles in the final state, there are two entries per 
event in each plot. 
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Fig. 4. Distributions of difference in log-likelihood from ensembles of simulated 
samples generated according to model B (left distribution) and model A (right 
distribution). 
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